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Abstract: Linear optimization problems are investigated whose parameters 
are uncertain. We apply coherent distortion risk measures to capture the pos- 
sible violation of a restriction. Each risk constraint induces an uncertainty 
set of coefficients, which is shown to be a weighted-mean trimmed region. 
Given an external sample of the coefficients, an uncertainty set is a convex 
polytope that can be exactly calculated. We construct an efficient geometri- 
cal algorithm to solve stochastic linear programs that have a single distortion 
risk constraint. The algorithm is available as an R-package. Also the algo- 
rithm's asymptotic behavior is investigated, when the sample is i.i.d. from 
a general probability distribution. Finally, we present some computational 
experience. 
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1 Introduction 



Uncertainty in the coefficients of a linear program is often handled by prob- 
ability constraints or, more general, bounds on a risk measure. The random 
restrictions are then captured by imposing risk constraints on their violation. 
Consider the linear program 



c'x 



mm 



s.t. Ax > b . 



(1) 



and assume that A is a stochastic m x d matrix and b e IR m . This is a 
stochastic linear optimization problem. To handle the stochastic side condi- 
tions a joint risk constraint, 

p m (Ax-b)<0, (2) 

may be introduced, where p m is an m-variate risk measure. E.g. with 
p m (Y) = Prob[y < 0] — a the restriction ^ becomes 

Prob[Ax > b] > 1 - a , (3) 

and a usual chance- constrained linear program is obtained. Alternatively, the 
side conditions may be subjected to separate risk constraints, 



bj)<0, j = l 



m 



(4) 



with Aj denoting the j-th row of A. In (|4j) each side condition is subject 
to the same bound that limits the risk of violating the condition. A linear 
program that minimizes c'x subject to one of the restrictions ^ or ^ is 
called a risk- constrained stochastic linear program. 

For stochastic linear programs (SLPs) in general and risk-constrained SLPs 



in particular, the reader is e.g. referred to Kail and Mayer (2010). What we 



call a risk measure here is mentioned there as a quality measure, and useful 
representations of the corresponding constraints are given. As most of the 



literature, Kail and Mayer (2010) focus on classes of SLPs with chance con- 



straints that lead to convex programming problems, since these have obvious 
computational advantages; see also Prekopa (1995). Our choice of the qual- 



ity measure, besides its generality, enjoys a meaningful interpretation and, 
as will be seen later, enables the use of convex structures in the problem. 
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In the case of a single constraint (m = 1) notate 

p(a'x - b) < . (5) 

A practically important example of an SLP with a single risk constraint 
^ is the portfolio selection problem. Let fx, .. . ,fd be the return rates on 
d assets and notate f = (fi, . . . , f^)'. A convex combination of the assets' 
returns is sought, Fx = Ylj=i^j x 3f that has maximum expectation under a 
risk constraint and an additional deterministic constraint, 

max £?[f]'x, s.t. p(r'x) < p , x G C , (6) 

xgC 

where p is a risk measure, po G 1R is a given upper bound of risk (a nonneg- 
ative monetary value), and C G IR d is a deterministic set which restricts the 
coefficients Xk in some way. For example, if short sales are excluded, C is 
the positive orthant in M. d . The solution x* is the optimal investment under 
the given model. We will see that, if a solution exists, it is regularly finite 
and unique. In our geometric approach such a solution corresponds to the 
intersection of some line and a convex body that both contain the point E[r]. 

Regarding the choice of p, two special cases are well known. First, let p(r'x) = 
Prob[r'x < — Vq] and po = a. Then the optimization problem (J6]) says: 
Maximize the mean return Z?[f'x] under the restrictions x G C and 

V@R a (f'x) < v . 

That is, the value at risk V@R a of the portfolio return must not exceed the 
bound v . Second, let 

p(f'x) = -- ! Qr>x(t)dt, (7) 
a Jo 

where Qz signifies the quantile function of a random variable Z. This means 
that the expected shortfall of the portfolio return is employed in the risk 
restriction. 

In practice, a has to be estimated from data. If the solution of the SLP is 
based on a sample of observed coefficient vectors a 1 , . . . , a n G M. d , that is, on 
an external sample, the SLP is mentioned as an empirical risk- constrained 
SLP. In other words, we assume that a follows an empirical distribution that 
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gives equal mass - to some observed points a 1 ,..., a" £ M. d . Rockafellar 



and Uryasev (2000) investigate an empirical stochastic program that arises in 



portfolio choice when the expected shortfall of a portfolio is minimized. They 
convert the objective into a function that is convex in the decision vector x 
and optimize it by standard methods. This approach is commonly used in 
more recent works of these and other authors on portfolio optimization. 



A more complex situation is investigated by Bertsimas and Brown (2009) 



who discuss the risk-constrained SLP with arbitrary coherent distortion risk 
measures, which also include expected shortfall. These allow for a sound 
interpretation in terms of expected utility with distorted probabilities. For 
the linear restriction a so called uncertainty set is constructed which consists 



of all coefficients satisfying the risk constraint. Bertsimas and Brown (2009) 



discuss the uncertainty set that turns the SLP into a minimax problem, called 
robust linear program; however they provide no optimal solution of this pro- 
gram. The uncertainty set is a convex body and, as it will be made precise 
below in this paper, comes out to equal a weighted-mean trimmed region. 



Natarajan et al. (2009), on the reverse, construct similar risk measures from 



given polyhedral and conic uncertainty sets. Pffug (2006) has proposed an 
iterative algorithm for optimizing a portfolio using distortion functionals, on 
each step adding a constraint to the problem and solving it by the simplex 
method. Meanwhile, many other authors have recently contributed to the 
development of robust linear programs related to risk-constrained optimiza- 



tion problems, see, e.g. Nemirovski and Shapiro (2006), Ben-Tal et al. (2009) 



and Chen et al. (2010) For a review of robust linear programs in portfolio 



optimization the reader is referred to Fabozzi et al. (2010). 



In this paper we contribute to this discussion in three respects: 



1. The uncertainty set of an SLP under a general coherent distortion risk 
constraint is shown to be a weighted-mean region, which provides a 
useful visual and computable characterization of the set. 

2. An algorithm is constructed that solves the minimax problem over the 
uncertainty set, hence the SLP. 

3. If the external sample is i.i.d. from a general probability distribution, 
the uncertainty set and the solution of the SLP are shown to be con- 
sistent estimators of the uncertainty set and the SLP solution. 
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The paper is organized as follows: In Section [2] constraints on distortion risk 
measures and their equivalence to uncertainty sets in the parameter space are 
discussed; further these uncertainty sets are shown to be so called weighted- 
mean trimmed regions that satisfy a coherency property. In Section [3] a 
robust linear program is investigated by which the SLP with a distortion risk 
constraint is solved. Section [4] introduces an algorithm for this program and 
discusses sensitivity issues of its solution. In Section [5] we address the SLP 
and its solution for generally distributed coefficients and investigate the limit 
behavior of our algorithm if based on an independent sample of coefficients. 
Section [6] contains first computational results and concludes. 

2 Distortion risk constraints and weighted- 
mean regions 

Let us consider a probability space {Q, J 7 , P) and a set 1Z of random variables 
(e.g. returns of portfolios). A function p : 7Z — > W is a law invariant risk 
measure if for Y, Z G 1Z it holds: 

1. Monotonicity: If Y is pointwise larger than Z then it has less risk, 
P(Y) < p(Z) . 

2. Translation invariance: p(Y + 7) = p{Y) — 7 for all 7 G R . 

3. Law invariance: If Y and Z have the same distribution, Py = Pz, then 
p{Y) = p(Z) . 

p is a coherent risk measure if it is, in addition, positive homogeneous and 
subadditive, 

4. Positive homogeneity: p{XY) = Xp(Y) for all A > , 

5. Subadditivity: p(Y + Z) < p(Y) + p(Z) for all Y, Z e TZ . 

The last two restrictions imply that diversification is encouraged - a crucial 
property for the risk management. Distortion risk measures are essentially 
the same as spectral risk measures. For the theory of such risk measures, 
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see e.g. Follmer and Schied (2004). A function p : 1Z — > K is said to satisfy 
the Fatou property if liminf n ^ l00 p(F n ) > p(Y) for any bounded sequence 
converging pointwise to Y . With the notion of coherent risk measures, we 



reformulate a fundamental representation result of Huber (1981): 



Proposition 1. p is a coherent risk measure satisfying the Fatou property if 
and only if there exists a family Q of probability measures that are dominated 
by P (i.e. P(S) = => Q(S) = for any S G J= and Q G Q) such that 

p(Y) = sup E Q (-Y) . 



We say that the family Q generates p. In particular, let (Q, A) = (M. d , B d ) 
and P be the probability distribution of a random vector a. Huber's Theorem 
implies that for any coherent risk measure p there exists a family G of P- 
dominated probabilities on B d so that 

p(a'x - b) < 



p(a'x) < -b 

inf E G (a'x) > b 

GeG 

E G (a'x) > b for all GeG. 



Let us denote by A n the unit simplex in IR n , 

n 

A n = {x e W 1 1 x k = h %k > VA;} . 
fe=i 

Then, if a has an empirical distribution on n given points in M. d , any subset 
Q of A n corresponds to a family of P-dominated probabilities, and thus 
defines a coherent risk measure p. As an immediate consequence of Huber's 
theorem an equivalent characterization of the risk constraint is obtained (see 
also Bertsimas and Brown (2009)): 

Proposition 2. Let p : 1Z — > K be a coherent risk measure and let a have an 
empirical distribution on a 1 , ... , a n e R d . Then there exists some Q p C A n 
such that 

p(a'x — 6) < a'x > b for all 

aeU p := conv{a G M. d \ a = [a 1 , . . . , a n ]q | q G Q p } . 
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Here, conv(W) denotes the convex closure of a set W . Proposition [2] says 
that a deterministic side condition a'x > b holding uniformly for all a in 
the uncertainty set IA P is equivalent to the above risk constraint (|5| on the 
stochastic side condition. This will be used below in providing an algorithmic 
solution of the risk-constrained SLP. 



2.1 Distortion risk measures 



A large and versatile subclass of risk measures is the class of distortion risk 



measures (Acerbi, 2002). Again, let Qy denote the quantile function of a 



random variable Y. 

Definition 1 (Distortion risk measure). Let r be an increasing function 
[0, 1] — > [0, 1]. The risk measure p given by 



P(Y) 



Q Y {t)dr{t) 



(8) 



is a distortion risk measure with weight generating function r. 



A distortion risk measure is coherent if and only if r is concave. For ex- 
ample, with r(t) — if t < a and r{t) — 1 if t > a, the value at risk 
V@K a (Y) = —Qy{ol) is obtained, which is a non-coherent distortion risk 
measure. A prominent example of a coherent distortion risk measure is the 
expected shortfall, which is yielded by r(t) — t/a if t < a and r(t) = 1 oth- 
erwise. Note that with r(t) = t, the risk measure becomes the expectation 
of —Y. A general distortion risk measure p(Y) can thus be interpreted as 
the expectation of — Y with respect to a probability distribution that has 
been distorted by the function r. In particular, a concave function r distorts 
the probabilities of lower outcomes of Y in positive direction (the lower the 
more) and conversely for higher outcomes (the higher the less). In empirical 
applications, coherent distortion risk measures other than expected shortfall 



have been recently used by many authors; see, e.g., Adam et al. (2008) for a 



comparison of various such measures in portfolio choice. 

An equivalent characterization of a coherent distortion risk measure is that 



it is coherent and comonotonic (Acerbi (2002)). p is comonotonic if 



p(Y + Z) = p(Y) + p(Z) for all Y and Z that are comonotonic, 
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i.e., that satisfy (Y(cj) - Y(u'))(Z(u) - Z(u')) > for every u,u' G Q. 
If Y has an empirical distribution on yi, . . . , y n £ R, the definition ^ of a 
distortion risk measure specializes to 

n 

P(>0 = -I>fN' (9) 
i=i 

where yui are the values ordered from above and qi are nonnegative weights 
adding up to 1. (Observe that qi = r(yr n+i-» i) — rfyrn-Q.) Then, the distor- 
tion risk measure pi) is coherent if and only if the weights are ordered, i.e. 
q g A™ := {q G A n | < q 1 < ■ ■ ■ < q n }. 



2.2 Weighted-mean regions as uncertainty sets 



If p is a coherent distortion risk measure, the uncertainty set U p has a special 
geometric structure, which will be explored now in order to visualize the 
optimization problem and to provide the basis for an algorithm. We will 
demonstrate that U p equals a so called weighted-mean (WM) region of the 
distribution of a. 

Given the probability distribution Fy of a random vector Y in IR d , weighted- 
mean regions form a nested family of convex compact sets, {D a (Fy)} a e[o,i]i 
that are affine equivariant (that is D a (FAx+b) — A D a (F Y ) + b for any regular 
matrix A and b G lR d ). By this, the regions describe the distribution with 
respect to its location, dispersion and shape. Weighted-mean regions have 



been introduced in Dyckerhoff and Mosler (2011) for empirical distributions 



and in Dyckerhoff and Mosler (2012) for general ones. 

For an empirical distribution on a 1 , . . . ,a n e IR d , a weighted-mean region is 
a polytope in IR d and defined as 



conv 



E 



t(3) 



7T permutation of {1, . . . , n} 



(10) 

Here w a = [w a> i, . . . , iu a>n ]' is a vector of ordered weights, i.e. w a G A<, 
indexed by < a < 1 that for a < (3 satisfies 



3=1 3=1 



J ■ 



VA; 



,n. 



(11) 



S 



Any such family of weight vectors {w Q }o< a <i specifies a particular notion 
of weighted-mean regions. There are many types of weighted-mean regions. 
They contain well known trimmed regions like the zonoid regions, the ex- 
pected convex hull regions and several others. For example, 



J_ 

not 


if 3 


> n — 


[na\ , 


net— [net] 
na 


if 3 


= n — 


\na\ , 





if 3 


< n — 


\na\ , 



< a < 1, defines the zonoid regions. However some popular types of 
trimmed regions, such as Mahalanobis or halfspace regions, are no weighted- 
mean regions. 

A WM region is characterized by its projections on lines. Note that each 
p G S^ 1 , where S^ -1 is a (d — l)-variate unit sphere, yields a projection of 
the data a 1 , . . . , a" on the line generated by p and thus induces a permutation 
7T P of the data, 

p' a ^p(i) < p V p(2) < • - - < p'a'pW . 

The permutation is not necessarily unique, and let //"(a 1 , . . . , a n ) denote the 
set of all directions p e that induce a non-unique permutation n p . 

Dyckerhoff and Mosler ( |2011 ) have shown that the support function h a of 



D Vfa = /^^(a 1 , . . . , a n ) amounts to 



Up) = w ^p'^ p(i) , p e S d -' . (12) 



It follows that, whenever ir p is unique, the polytope D Wa has an extremal 
point in direction p, which is given by 

n n 

£ WajaT*® = £ "a,*?®** - P e S"- 1 \ #(a\ . . . , a") . (13) 
j=i i=i 



Now we are moving to the main result of this section, which will be the 
Theorem |TJ From (|9j) and (jl3j) it is seen that, with % = w aiT -\^ and = 
— p'a\ the extreme point of the projection of D Wa on the p-line is obtained 
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by applying a q-distortion risk measure to the projected data points. Now, 
setting 

Q P = {q e A" I q = K j7r - 1(1) , . . . w a ^ (n) ), p G S^ 1 \ H{*\ . . . , a")} 

and Up = conv{a | a = Y^i=\ Qi^i ^ ^ 2p) we obtain that all extreme points 
of D Wa are in U p , hence D Wq C U p . On the other hand, for every q G Q p 
it holds that Y17=i e ^ w <* > which implies W p C _D Wq . We conclude 
U p = D Wa . ^ 

Thus we have proven the equality between the distortion risk constraint 
feasible set and a properly chosen WM region, which is a reformulation of 



Theorem 4.3 from Bertsimas and Brown (2009): 
Theorem 1. 

{x G R d | p(a'x -6)<0} = {xGl d |a'x>&VaG Avja 1 , . . . , a n )}, (14) 
where a 1 , . . . , a n is an external sample of the parameter vector a. 

Recall that D Wa (a 1 , . . . , a n ) is a <i-dimensional convex polytope, and thus 
the convex hull of a finite number of points (its vertices) or, equivalently, 
a bounded nonempty intersection of a finite number of closed halfspaces 
(that contain its facets). By this the calculation and representation of such a 
polytope can be done in two ways: either by its vertices or by its facets. Recall 
that a nonempty intersection of the polytope's boundary with a hyperplane 
is a facet if it has an affine dimension d — 1, and a ridge if it has an affine 
dimension d — 2. It is called an edge if it is a line segment, and a vertex 
if it is a single point. In general, each facet of a polytope in M. is itself a 
polytope of dimension d — 1 and has at least d vertices. With WM regions 
the number of a facet's vertices can vary considerably; it ranges between d 



and d\ (Bazovkin and Mosler, 2012). That is why in calculating WM regions 



a representation by facets is preferable. 

The vertices of a polytope are its extreme points. From above we know that 
the directions p G S d ~ 1 \H(sL 1 , . . . , a n ) belong to vertices, while the directions 
p G H^sl 1 , . . . , a n ) belong to parts of the boundary that have affine dimension 
> 1. 

In the context of risk measurement it is crucial that the WM regions pos- 
sess two properties that enable them to generate coherent risk measures: 
monotonicity and subadditivity. 
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Proposition 3. (Coherency properties of WM regions) 



1. Monotonicity: Ifzk < yk holds for all k (in the componentwise ordering 
ofR d ), then 

Av Q (yi, • • • ,y n ) C D Wa (zi,...,z n ) ©IR+, and 
D Wa ( Zl ,...,z n ) C D Wa ( yi ,...,y n )®R d _. 

2. Subadditivity: 
Av«(yi + zi, . . . ,y n + z„) c f Wa (yi, . . . ,y„) © Av Q (zi, . . . ,z n ) . 

In this Proposition the symbol © is the Minkowski addition, A(& B = {a + 



b\a € A, b e -B} for A and B C M d . For a proof, see Dyckerhoff and Mosler 



(2011). 



The subadditivity property of WM regions is an immediate extension of 
the subadditivity restriction usually imposed on univariate risk measures. 
In dimensions two and more it has an interpretation as a dilation of one 
trimmed region by the other. To understand this better let us consider the 
simple example of Minkowski addition given in Figure [T] The figure exhibits 
a solid triangle with one vertex at the origin and a dotted-border quadrangle. 
Now move the triangle in a way that its lower left corner passes all points 
of the quadrangle. At each point of the quadrangle we get a copy of the 
initial triangle (with a dashed border) shifted by coordinate of the point. 
The union of all these triangles gives us the Minkowski sum of the initial two 
sets, which is the big heptagon in the picture. Observe that, if the rectangle 
is moved around the triangle, the same sum is obtained. The subadditivity 
states that if, e.g., these two figures are WM regions -D Wa (yi, • . . , y n ) and 
L> Wq (zx, . . . , z n ) respectively, the -D Wa (yi + zi, . . . , y n + z n ) is contained by 
the heptagon. 
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Figure 1: An illustration of the subadditivity property. 



3 Solving the SLP with distortion risk con- 
straint 

3.1 Calculating the uncertainty set 

In the previous section we have shown that the uncertainty set U p equals the 
weighted- mean (WM) region D Wa for a properly chosen weight vector w a . 
Bazovkin and Mosler| (2012) provide an algorithm by which this WM region 



can be exactly calculated in any dimension d. The results can be visualized 
in dimensions two and three; for examples see Figure [2] 



It has been also shown in Bazovkin and Mosler (2012) that the number of 



vertices of a facet can be as much as d\ . Therefore the representation of a 
WM region by its vertices appears to be less efficient than that by its facets. 
In the sequel, we will use the facet representation for solving the SLP. 



3.2 The robust linear program 

The robust linear program to be solved is 

c'x — > min s.t. a'x > b for all a £ U , (15) 
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Figure 2: Visualization of WM regions by the R package WMTregions. Left 
panel: Facets of a three-dimensional region in M. 3 . Right panel: Vertices of a 
four-dimensional region projected on a subspace of R 3 . 



where the subscript p has been dropped for convenience. The side condition 
is rewritten as 



Note that X, as a weighted-mean region, is a convex polyhedron. So, a 
linear goal function is to be minimized on a convex polyhedron. Obviously, 
any optimal solution will lie on the surface of X. 

3.3 Finding the optimum on the uncertainty set 

In constructing an algorithm for the robust linear program, we will explore 
the set X of feasible solutions and relate it to the uncertainty set U in the 
parameter space. It will come out that the space of solutions x and the space 
of coefficients a are, in some sense, dual to each other. The following two 
lemmas provide the connection between X and U. First we demonstrate that 
X is the intersection of those halfspaces whose normals are extreme points 




X a = {x|a'x > b} . 



(16) 



of U. 
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Lemma 1. It holds that 

X = P| {x I a'x > b} = p| {x I a'x > b} . 

Proof. We show that flasextw^a c X " for all ufW; then Daeextw^a C 
f"| agW X a . The opposite inclusion is obvious. Assume ufW. Then, as U is 
convex and compact, u is a convex combination of some points a 1 , . . . , a £ G 
extU, i.e. u = Ylt=i ^« ai "with Aj > and ^1=1^ = 1' anc ^ f° r an y x e 
Has ext w holds x G X a3 and a J 'x > b for all j, hence u'x = Ylt=i Aja j/ x > 
b, that is, x G X u . □ 

Lemma [T] says that each facet of the set X of feasible solutions corresponds 
to a vertex of the uncertainty set U. Hence it is sufficient to consider the 
extreme points of the uncertainty set. 

As a generalization of Lemma [TJ we may prove by recursion on k: Each 
k- dimensional face of the feasible set corresponds to a (d — A;)-dimensional 
face of the uncertainty set in the solution space. This resembles the dual 



correspondence between convex sets and their polars (cf. e.g. Rockafellar 



( 1997| )). However, in contrast to polars, in our case the correspondence is 



not reflexive. 

From Lemma [T] it is immediately seen, how the robust optimization problem 
contrasts with a deterministic problem, where the empirical distribution of 
a concentrates at some u G U. Observe that the deterministic feasible set 
X is just a halfspace, X u = {x | u'x > b}. In the general robust case a 
halfspace is obtained for each a G extU, and the robust feasible set X is 
their intersection. The halfspaces are bounded by hyperplanes with normals 
equal to a G extW, and their intercepts are all the same and equal to b. 
Consequently, the robust feasible set X is always included in the deterministic 
feasible set X u , 

X C X u for any uGW. 

Moreover, the two feasible sets cannot be equal unless each element of U is a 
scalar multiple of u with a factor greater than one, hi C {a | a = Au, A > 1}. 
Consequently, the minimum value of the robust stochastic LP cannot be 
smaller than the value of an LP with any deterministic parameter u chosen 
from the uncertainty set. Figure |3] (left panel) illustrates how a deterministic 
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deterministic Xa 



x-space a-space 
(Solution space) (Parameter space) 

Figure 3: Deterministic and robust cases: feasible set (left panel), uncertainty 
set (right panel). 

feasible set in dimension two compares to a general robust one: The line that 
bounds the halfspace X u 'folds' into a piecewise linear curve delimiting X . 

Let 

U x = {a G R d |a' x > b} , x G R d . 
Lemma 2. It holds that 

x<=X xe ext X 

Moreover, each vertex x G ext X corresponds to a facet oflA. 

Proof. By Lemma [T] we have x G X ^ a'x > b for all aGW. Now let 
aeW; then for any x G X it holds that a'x > b, hence a G £/ x . Conclude 
U. C fixe A" ^x- Further, it is clear that an extreme point x G ext X yields a 
facet of IA. □ 

Remark. While U is always compact, X is in general not. Therefore neither 
inclusion holds with equality. 

The ordinary simplex algorithm, operating on the vertices of X, constructs 
a chain of adjacent facets in the space of parameters. The chain ends at 
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Solution space 



Parameter space 



Figure 4: Duality between spaces. 



the solution of the optimization task. Notice that this chain corresponds 
to a chain of facets of the uncertainty set. So, in principle we could try to 
calculate this chain of facets in the parameter set. However, in our algorithm, 
another way is pursued to find the optimal solution. 

To manage this task let us consider the goal function c'x. In the parameter 
space c corresponds to a point or a direction. In the solution space it corre- 
sponds to all hyperplanes that have c as their normal. To produce all these 
hyperplanes in the parameter space, c has to be multiplied with some scaling 
factor. Hence the hyperplanes are obtained by passing through a straight 
ray ip starting at the origin and containing c. 

Next we search the intersection of U with the ray (p. Note that finding 
the intersection of a line and a polyhedron in IR 3 is an important problem 



in computer graphics (cf. Kay and Kajiya (1986)). The same principle is 



employed for a general dimension d. The uncertainty set U is the finite 
intersection of halfspaces Hj, j = 1 . . . J, each being defined by a hyperplane 
Hj with normal n., pointing into Hj and an intercept dj. 

Consider some point u on the ray ip that is not in U. Compute for all 
halfspaces Hj that do not include u, i.e. where (u'n,- — dj) < holds. (In 
other words, Hj is visible from u.) Find j* at which this value is the largest. 
Recall that moving a point u along ip is equivalent to multiplying u by some 
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constant. The furthest move is given by the biggest constant. The optimal 
solution x* of the robust SLP has to satisfy a'x* > b, which is equivalent to 



Besides the regular situation described above, two special cases can arise: 

1. There is no facet visible from the origin. This means that no solution 
is obtained. 

2. tp does not intersect U. Then the whole procedure is repeated with the 
opposite ray —tp. If this still gives no intersection, an infinite solution 



Finally, we like to point out that not the whole polytope U needs to be 
calculated but only a part of it which intersects the ray tp. In searching 
for the optimum not all F facets need to be checked, but only a subset of 
the surface where the intersection will happen. Such a filtration makes the 
procedure more efficient. The search for a proper subset can be driven by 
geometrical considerations. Let x* be an optimal solution of the robust SLP. 
A subset U c s of U will be mentioned as an efficient parameter set if 

• x * e rw ff { x i a/x ^ h > c x and 

• a, d e U e s, a'x > b, d'x > b implies a = d . 

That is to say, U e s is the minimal subset of U containing all facets that can 
be optimal for some c. 

Proposition 4. U e g is the union of all facets ofU for which dj > holds. 




Hence, to obtain x*, the normal has to be scaled with the constant ^ 




(17) 



exists. 
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In other words, an efficient parameter set U e s consists of that part of the 
surface of U that is visible from the origin 0. The proof is obvious. 

To visualize the efficient parameter set we will use the augmented uncertainty 
set, which is defined as 

{a|a = Aa*, A > 1, a* e U eS } . 

It includes all parameters that are dominated by U c s] see the shaded area in 
the right panel of Figure [5] 

So far we have assumed that b > 0. It is easy to show, that with b < 
we have to construct the intersection of ip with the part of the surface of hi 
that is invisible from the origin 0, which is U e s in this case. In the sense of 
Proposition EJ U e g contains all facets of U with dj < 0. Obviously, U e s is 
always non-empty in this case, which, in turn, means that the existence of 
a solution is guaranteed. However, the solution can be infinite if ip does not 
intersect U e s. 

The situation of b < is common in the maximizing SLPs. Really, if we have 
the model 

c'x — > max s.t. a'x < b for all a G U , (18) 
it is possible to rewrite it as follows: 

(-c)'x — > min s.t. (-a)'x > -b for all a e U . (19) 



Clearly, (19) is equivalent to (15) except of the negativity of the coefficient b. 



4 The algorithm 

In this part an accurate procedure of obtaining the optimal solution is given. 
Input: 

• a vector c 6 M. d of coefficients of the goal function, 

• an external sample {a 1 , . . . ,a n } C M. d of coefficient vectors of the re- 
striction, 
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a right-hand side b G K of the restriction, 

a distortion risk measure p (defined either by name or by a weight 
vector). 



Output: 

• the uncertainty set U of parameters given by 

— facets (i.e. normals and intercepts), 

— vertices, 

• the optimal solution x* of the robust LP and its value c'x*. 
Steps: 

A. Calculate the subset ZY e fr C U consisting of facets {(n?, dj)}j e j. 

B. Create a line p passing through the origin and c. 

C. Search for a facet Hj t of U c s that is intersected by (p: 

a. Select a subset U se i C U eS of facets: This may be either U e s itself 
or its part where the intersection is expected; U se i = {(n^, dj) \ j G 
J se i}. For example, we can search the best solution on a pre-given 
subset of parameters. The other possible filtration is iterative 
transition to a facet with better criterion value. 

b. Take a point u = Ac, A > 0, outside the augmented uncertainty 
set. Find the j* = argmax{A 3 - = ^r|Aj > 0} jeJsei cj- For the 

3 3 

case b < just replace argmax with argmin. 
I. If ip does not intersect U e g, then the solution is infinite. If 



b > 0, then repeat |C.b.| for the opposite ray — else stop. 

II. If in the case b > the inner part of U contains the origin, 
then no solution exists; stop. 

c. x* = ^-n^ is the optimal solution of the robust LP. 
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Solution space 



Parameter space 



Figure 5: Finding the optimal solution on the uncertainty set. 



In fact, the line (p consists of points that correspond to hyperplanes whose 
normal is the vector c in the dual space. One part of ip is dominated by 
points from U e g, while the other is not (which results from Proposition [1]). 
The crossing point a* defines the hyperplane that touches the feasible set at 
the optimum as its dual. 

Moreover, a typical nonnegativity side constraint x > can be easily ac- 
counted for in the algorithm. In considering this, the search for facets has 
just to be restricted to those having nonnegative normals. 

To solve the portfolio selection problem ^ with the algorithm, we treat the 
realizations of the vector — f of losses rates as {a 1 , . . . ,a n }, and minimize 
c'x with c = - Yli=i a • This corresponds to transforming the maximizing 



SLP by (19) and running the above outlined procedure. Note that both tp 
and U contain the point - Y^Ji=i a *' that ls i they always intersect, which, in 
turn, guarantees the existence of a finite solution. To meet a unit budget 
constraint, the solution x* is finally scaled down by Ylj=i x j — 1- Recall that 
the risk measure is, by definition, scale equivariant. 
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4.1 Sensitivity and complexity issues 



Next we like to discuss how the robust SLP and its optimal solution behave 



when the data {a 1 , . . . , a n } on the coefficients are slightly changed. From ( 12 ) 
it is immediately seen that the support function hu of the uncertainty set is 
continuous in the data aP as well as in the weight vector w Q . (Note that the 
support function hu is even uniformly continuous in a 1 , . . . , a™ and w a , which 
is tantamount saying that the uncertainty set U is Hausdorff continuous in 
the data and the risk weights.) Consequently, a slight perturbation of the 
data will only slightly change the value of the support function of U, which 
is a practically useful result regarding the sensitivity of the uncertainty set 
with respect to the data. The same is true for a small change in the weights 
of the risk measure. 

We conclude that the point a J * where the line through the origin and c cuts 
U depends continuously on the data and the weights. However this is not 
true for the optimal solution x*, which may 'jump' when the cutting point 
moves from one facet of U to a neighboring one. 

The theoretical complexity in time of finding the solution is compounded 
from the complexity of one transition to the next facet and by the whole 



number of such transitions until the sought-for facet is achieved. Bazovkin 



and Mosler (2012 ) have shown that the transition has a complexity of 0(d 2 n). 
In turn, in the same paper the number of facets N(n, d) of an WM region is 
shown to lie between 0(n d ) and 0(n 2d ) depending on the type of the WM 
region. Thus, it is easily seen, that an average number of facets in a facets 
chain of a fixed length is defined by the density of facets on the region's 
surface, y/N(n, d), and is estimated by a function between 0{n) and 0{n 2 ). 
The overall complexity is then 0(d 2 n 2 ) up to 0(d 2 n 3 ). Notice, that the 
lower complexity is achieved for zonoid regions, namely when the expected 
shortfall is used for the risk measure. 



4.2 Ordered sensitivity analysis 

Also alternative uncertainty sets may be compared that are ordered by inclu- 
sion. From Lemma [T] it is clear that the respective sets of feasible solutions 
are then ordered in the reverse direction; see e.g. Figure |6j In particular 
we may consider the robust LP for two alternative distortion risk measures 
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Reversed central regions 
forx 



Central regions 
for a 



Figure 6: Example of the 'reversed' central regions in the dimension 2. 



which are based on weight vectors w Q and wp, respectively, that satisfy the 
monotonicity restriction ( 11 ). Then the resulting uncertainty sets are nested, 
Up C U a and so are, reversely, the feasible sets, Xp D X a . This is a useful 
approach for visualizing the sensitivity of the robust LP against changes in 
risk evaluation. 



5 Robust SLP for generally distributed coef- 
ficients 

So far an SLP (II]) has been considered where the coefficient vector a follows 
an empirical distribution. It has been solved on the basis of an external sam- 
ple {a 1 , . . . , a n }. In this section the SLP will be addressed with a general 
probability distribution P of a. We formulate the robust SLP in the gen- 
eral case and demonstrate that the solution of this SLP can be consistently 
estimated by random sampling from P. 

Consider a distortion risk measure p pi) that measures the risk of a general 
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random variable Y and has weight generating function r, p(Y) = — Qy(t)dr{t). 
Similarly as in Section 2.2 a convex compact U in W 1 is constructed through 
its support function hy, 



MP J 



Q p/ ~ a {t)dr{t) . 



Now, let a sequence (a n ) neN of independent random vectors be given that 
are identically distributed with P, and consider the sequence of random un- 
certainty sets U n based on a 1 , . . . , a 
shown: 



Dyckerhoff and Mosler (2011) have 



Proposition 5 (Dyckerhoff and Mosler (2011)). U n converges to U almost 
surely in the Hausdorff sense. 



The proposition implies that by drawing an independent sample of a and 
solving the robust LP based on the observed empirical distribution a con- 
sistent estimate of the uncertainty set U is obtained. Moreover, the cutting 
point a- 7 *, where the line through the origin and c hits the uncertainty set, is 
consistently estimated by our algorithm. But, in particular for a discretely 
distributed a, the optimal solution x* need not be a consistent estimate, as 
it may perform a jump when aJ* moves from one facet of hi to neighboring 
one. 



6 Concluding remarks 



A stochastic linear program (SLP) has been investigated, where the coeffi- 
cients of the linear restrictions are random. Risk constraints are imposed 
on the random side conditions and an equivalent robust SLP is modeled, 
whose worst-case solution is searched over an uncertainty set of coefficients. 
If the risk is measured by a general coherent distortion risk measure, the 
uncertainty set of a side condition has been shown to be a weighted-mean 
region. This provides a comprehensive visual and computable characteriza- 
tion of the uncertainty set. An algorithm has been developed that solves the 
robust SLP under a single stochastic constraint, given an external sample. It 
is available as an R-package StochaTR (Bazovkin and Mosler, 2011). More- 
over, if the data are generated by an infinite i.i.d. sample, the limit behavior 
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of the solution has been investigated. The algorithm allows the introduc- 
tion of additional deterministic constraints, in particular, those regarding 
nonnegativity. 



d\n 


1000 


2000 


3000 


4000 


5000 


10000 


15000 


20000 


25000 


3 


0.3 


1.14 


1.76 


2.92 


3.41 


6.18 


12.61 


15.06 


47.54 


4 


0.66 


2.21 


3.47 


4.48 


4.27 


7.68 


16.97 


20.04 




5 


1.85 


3.09 


5.68 


9.28 


11.03 


13.52 


27.34 


54.86 




6 


2.08 


4.41 


5.62 


14.99 


18.73 


25.07 


46.88 






7 


2.16 


6.22 


13.3 


25.44 


28.56 


52.33 








8 


4.18 


9.78 


20.18 


31.82 


34.23 










9 


5.18 


14.75 


24.11 


35.94 


61.14 










10 


6.17 


16.97 


33.82 


42.11 


67.06 











Table 1: Running times of StochaTR for different n and d (in seconds). 

Table [I] reports simulated running times (in seconds) of the R-package for the 
5%-level expected shortfall and different d and n. The data are simulated 
by mixing the uniform distribution on a <i-dimensional parallelogram with 
a multivariate Gaussian distribution. In light of the table the complexity 
seems to grow with d and n slower than 0(d 2 n 2 ). 

Besides this, we contrast our new procedure with the seminal approach of 



Rockafellar and Uryasev (2000), who solve the portfolio problem by optimiz- 
ing the expected shortfall with a simplex-based method. In illustrating their 
method, they simulate three-dimensional normal returns having specified ex- 
pectations and covariance matrices. We have applied our package to likewise 
simulated data. The results are exhibited in Table [2] For a comparison, some 



cells contain also a second value, which corresponds to the Rockafellar and 



Uryasev (2000) procedure and is taken from Table 5 there. 



a\n 


1000 


5000 


10000 


15000 


20000 


25000 


0.10 


1.1 (<5) 


7.2 (6) 


23.7 (20) 


46 


56.3 (45) 


74.4 


0.05 


0.5 (<5) 


4.7 (6) 


14.0 (12) 


20.0 


39.8 (40) 


53.2 


0.01 


0.3 (<5) 


2.3 (6) 


3.8 (6) 


7.9 


22.1 (50) 


38.5 



Table 2: Running times of StochaTR for different n and a (in seconds); in 



parentheses running times of Rockafellar and Uryasev (2000) 
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As we see from Table [2j the computational times of the two approaches do not 
much differ. However, our algorithm usually needs some dozens of iterations 



only, which is substantially less than the algorithm of Rockafellar and Uryasev 



( 2000[ ). Also, in contrast to the latter, where the resulting portfolio can vary 



between (0.42,0.13,0.45) for n = 1000 and (0.64,0.04,0.32) for n = 5000, 
we get a stable optimal portfolio. Our solution averages at (0.36,0.15,0.49), 
which has approximately the same V@R and expected shortfall as that in 
the compared study but yields a better value of the expected return. 

Finally, our approach turns out to be very flexible. In particular, non-sample 
information can be introduced into the procedure in an interactive way by ex- 
plicitly changing and modifying the uncertainty set. More research is needed 
in extending the algorithm to solve SLPs with multiple constraints (j2|. Also 
procedures that allow for a stochastic right-hand side in the constraints and 
a random coefficients in the goal function have still to be explored. 
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